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We develop methods for calculating the zero-frequency noise for quantum shuttles, i.e. nanoelec- 
tromechanical devices where the mechanical motion is quantized. As a model system we consider a 
three-dot array, where the internal electronic coherence both complicates and enriches the physics. 
Two different formulations are presented: (i) quantum regression theorem, and (ii) the counting 
variable approach. It is demonstrated, both analytically and numerically, that the two formulations 
yield identical results, when the conditions of their respective applicability are fulfilled. We describe 
the results of extensive numerical calculations for current and current noise (Fano factor), based on 
a solution of a Markovian generalized master equation. The results for the current and noise are 
further analyzed in terms of Wigner functions, which help to distinguish different transport regimes 
(in particular, shuttling vs. cotunneling) . In the case of weak inter-dot coupling, the electron trans- 
port proceeds via sequential tunneling between neighboring dots. A simple rate equation with the 
rates calculated analytically from the P(i5)-theory is developed and shown to agree with the full 
numerics. 

PACS numbers: 85.85.-|-j, 72.70.-|-m, 73.23.Hk, 73.63.-b 



I. INTRODUCTION 

As the advances of the technology push the size of the 
electronic components towards the atomic scale new in- 
teresting phenomena influencing the electronic transport 
emerge. New research fields, e.g. molecular electronics, 
spintronics, or nanoelectromechanical systems (NEMS) 
have appeared. A common theme is the combination of 
quantum transport and a subtle interplay between vari- 
ous degrees of freedom which plays an essential role for 
the functionality of the device. This paper focuses on 
the NEMSjii^ a logical extension of the now established 
technology of MEMS, where the electronic (or magnetic) 
degrees of freedom are coupled to a mechanical degree 
of freedom. While still in its infancy, NEMS have al- 
ready attracted much attention both experimentally^T^ 
and theoretically>i2i"i^ 

A measurement of the stationary IV-characteristic of 
a NEMS device does not always yield enough infor- 
mation to uniquely identify the underlying microscopic 
charge transport mechanism. A point in case is the 
Ceo SET experiment by Park et ali^ where two alterna- 
tive interpretations, namely incoherent phonon assisted 
tunnelingiSiSiiSSiSi or shuttling, are plausible. The 
current noise provides another important characteristics, 
supplementary to the mean current jSirJSi The Fano fac- 
tor, being the ratio between the zero-frequency compo- 
nent of the noise spectrum and the mean current, char- 
acterizes the degree of correlation between charge trans- 
port events and is a powerful diagnostic tool which helps 
to distinguish various transport mechanisms possibly re- 
sulting in the same mean current. Therefore, studies of 
the current noise in NEMS have formed an active field 
of researchiSSiSSiSiiMiS^ These studies considered noise in 



movable singe-electron transistors in a number of differ- 
ent configurations. 

To the best of our knowledge, the effects of internal co- 
herence of the electronic subsystem on the noise in NEMS 
have not been considered so far. The coherence is not a 
dominating feature in a system consisting of a single-level 
molecule or quantum dot. However, in a setup consisting 
of an array of dots the role of the electronic coherence 
within the array is of central importance. Its influence 
on the current in static quantum dot arrays has been 
studied intensivelji^^ and, more recently, also on the 
noise. Also, the mean current dependence on various 
system parameters in movable quantum dot arrays has 
already been studied>iSi2S Thus, the study of noise in a 
movable quantum dot array is the central theme in this 
work. 

Specifically, we study an array of three quantum dots in 
the strong Coulomb blockade regime with a movable cen- 
tral dot. This model was proposed as a quantum shuttle 
by Armour and MacKinnoni^ extending the original one- 
dot shuttle proposal by Gorelik et &\}^ The electronic 
coherence within the array combined with the mechani- 
cal degree of freedom changes qualitatively the transport 
through the array as compared to both a static array or 
a one-dot SET-NEMS. In particular, there are two com- 
peting electron transfer mechanisms through the array: 
either sequential tunneling or cotunneling (virtual tran- 
sition) via the central dot. The state of the oscillator fur- 
ther influences these two basic mechanisms which leads to 
a possibility of many different transport regimes depend- 
ing sensitively on the interplay of the parameters of the 
model. Roughly speaking, as we shall see cotunneling is 
associated with super-Poissonian values of the Fano fac- 
tor (sometimes as high as 50) while the sequential tun- 
neling is accompanied by sub-Poissonian Fano factors^ 



Similar conclusions have been reported in recent liter- 
ature for different but related systems, and a detailed 
discussion is given in sections to follow. 

We have recently published two Letters on quantum 
shuttleS(22iM and while the present paper addresses a 
somewhat different physical system, it makes heavy use 
of the techniques developed in the two Letters. Since we 
believe that the techniques may have a wide range of ap- 
plications, we use this opportunity to describe our general 
approach to quantum shuttles and expose the theoreti- 
cal machinery in more detail. The paper is organized 
as follows. In Section 2 we introduce our model of the 
three-dot quantum shuttle which is quite similar to the 
one considered in Ref. I (i. The total Hamiltonian con- 
sisting of the "system" (both mechanical and electronic 
degrees of freedom of the quantum dot array), the leads, 
and a generic heat bath is used to illustrate the derivation 
of a description based on Markovian generalized master 
equation which was the starting point of Ref. 16. Along 
the way from the Hamiltonian to the generalized mas- 
ter equation we identify several tacit assumptions used 
in previous studies (including ours) and point out sev- 
eral issues of potential importance not addressed so far 
within the field of NEMS. While we are not able to resolve 
all of these issues we believe that spelling them out is an 
important first step towards their solution. In particular, 
we address the problem of the assumed additivity of two 
kinds of baths acting on the system (the Fermi seas of 
the leads and the heat bath weakly coupled to the sys- 
tem). Another point of concern is the possible spurious 
breaking of the charge conservation by the weak-coupling 
prescription between the heat bath and the system with 
internal coherence. We close Sec. 2 with a short intro- 
duction to the superoperator formalism. 

In Section 3 we develop the theory of the zero- 
frequency component of the current noise spectrum for 
a NEMS device where the electron transfer between the 
system and the leads is described by a classical Markov 
process, i.e. in the wide band approximation and high 
bias limit. We present two methods of the evaluation of 
the noise spectra. If the whole system dynamics can be 
described by a Markovian generalized master equation 
we can use the quantum regression theorem. The other 
method relies on the counting variable approach and cal- 
culates the zero-frequency current noise as the charge 
diffusion coefficient across a given junction between the 
system and a lead. As we show further in Sec. 3 the two 
approaches yield equivalent results provided that the dy- 
namics of the system is (quantum) Markovian and that 
charge conserving approximations are used. We finish 
Sec. 3 by a qualitative discussion of the numerical evalu- 
ation of the noise spectra. This is a non-trivial task due 
to large dimensions of the involved matrices. Further de- 
tails of the numerical algorithm (Arnoldi iteration and 
generalized minimum residual method) are given in Ap- 
pendix A. 

We present the results of our numerical and analyti- 
cal calculations in Section 4. Generic features observed 
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FIG. 1: Schematic picture of the three dot system. The outer 
dots are fixed — the left one (L) at the position —xq and the 
right one (R) at xo, while the central one (C) can move (po- 
sition x) in a harmonic confining potential. It also interacts 
with a heat bath causing damping and thermal noise. The 
outer dots whose respective energy levels are de-aligned by 
the device bias (ei,) are coupled to the full or empty electronic 
reservoirs (leads), respectively. The current flows within the 
system due to tunneling between the left and central dot and 
the central and right dot and is described by the correspond- 
ing current operators Icl, Irc- 

in the numerical curves are interpreted phenomenologi- 
cally. Next, wc study different limiting cases. The first 
limit is that of small damping which is relevant for shut- 
tling accompanied by relatively small Fano factors (down 
to w 0.25) and strong inelastic cotunneling accompanied 
by huge Fano factors. These two mechanisms may coexist 
leading to a dramatic dependence of the Fano factor on 
parameters as the relative weight of the two mechanisms 
is changed. The second limit considered is the limit of 
weak coupling between adjacent dots which leads to se- 
quential tunneling assisted by an equilibrated oscillator, 
at least in a certain range of other parameters. In the se- 
quential tunneling limit we fully reproduce the numerical 
results with (semi-)analytic rate-equation-based theory 
with the rates determined by the standard P(i?)-theory 
as functions of the model parameters. The technical de- 
tails of the analytic calculations are sketched in Appendix 
B. We state our conclusions in Sec. 5. 



II. THREE-DOT QUANTUM DOT ARRAY 

A. The Model 

Armour and MacKinnon^^ introduced a model of a 
three-dot array whose central dot is movable. The ar- 
ray is assumed to be in the strong Coulomb blockade 
regime in which only two charge states (none or one ex- 
tra electron which we refer to as unoccupied or singly 
occupied) of the whole array, separated by an energy dif- 
ference Eo, are allowed in the considered bias range. This 
can be achieved by a suitable gating of the array which 
makes the two charge states energetically close while a 
very high charging energy prohibits addition or removal 
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of other electrons to/from the array. The array is coupled 
to two leads with a high bias applied between them. The 
bias is smaller than the charging energy for addition or 
removal of other electrons but otherwise it is the largest 
energy scale in the model. 

The moving central dot interacts with its surroundings 
and the dissipative dynamics is described by the interac- 
tion with a generic heat bath. We modify the original 
model slightly in that we do not consider the additional 
hard wall potential at the position of the outer dots ±xo 
employed by Armour and MacKinnon^^ so that the cen- 
tral dot moves in a strictly harmonic potential in our case 
(see Fig. ^ . While the hard wall potential is physically 
well motivated it complicates the numerical treatment 
and we believe that it does not have any significant im- 
pact on the nature of our results. Therefore, in our model 
the amplitude of oscillations in some regimes can exceed 
xo- The hard wall potential can be straightforwardly 
incorporated in our formalism. It should be noted, how- 
ever, that the various models for dissipation used in the 
literature, and also adopted in our work, are best justi- 
fied for the pure harmonic potential. Also, as in Ref.lll 
we consider spinless electrons. 

The Hamiltonian reads 



H — Hose + Hci + i/ol-osc + fflcads + ^^cl-lcads 



^^bath + H, 



bath -t- -fJ osc -bath + ^^CT 



where 



Ho 



IP 



2m 



(1) 



(la) 



describes the mechanical center-of-mass motion of the 
central dot as a one-dimensional harmonic oscillator with 
mass m and frequency ujq. The next two terms specify the 
electronic structure of the array in the strong Coulomb 
blockade regime (i.e. no double occupancy in the whole 
array — the vectors |/) with / = 0, L,C,R span its en- 
tire electronic Hilbert space) and the electromechanical 
coupling within the array 



Hci + Hci-c 



|L)(L|-^|i?)(i?|+eo 



+ tL{x) [\L){C\ + \C){L\ ) + tnix) {\C){R\ + \R){C\ 
-^^|C)(q 

(lb) 

with tL{x) = -Foe-"(^o+*), tR{x) = -Vbe"(*-^o). We 
associate the energies and eo with the left and 

right dot and the empty array, respectively, while the 
energy level of the central dot is chosen as the refer- 
ence energy, and hence put to zero. The device bias 
is the difference between the energy of the left and the 
right dot (which can be induced by suitable gating of the 
different dots) and 2xo is the distance between the two 
outer dots. The terms proportional to tL,R{x) describe 



a position-dependent hopping between the left and cen- 
tral or central and right dots enabling the tunneling cur- 
rent to flow through the array. These terms contribute 
both to the static part of the Hamiltonian (zeroth order 
in x) as well as to the electromechanical coupling. The 
parameter a equals the inverse tunneling length and de- 
termines the strength of the exponential ^-dependence 
of the hopping elements which may lead to the shuttling 
instabilityiifliiSiSS The last term gives the electromechan- 
ical coupling due to the electrostatic force acting on the 
oscillator when the central dot is charged. 

The outer dots of the array are assumed to couple via 
standard tunneling terms to two non-interacting leads: 



-fflcads + ^^ol-lcads — ^ £k[jc\pCkl3 
k-p=L,R 

+ E ykp(cip\m\ + \m\ckf3 



(Ic) 



The leads are held at different electrochemical poten- 
tials ^L.R whose difference gives the bias across the ar- 
ray. We assume that the tunneling densities of states 
r/3(e) = X Sfc l^fe/3p<5(e - £k(3) are energy-independent 
(and equal, just for convenience), i.e. r^(e) = F, known 
as the wide-band limit. It is necessary for the so-called 
first Markov approximation^^i^ used later on, to hold. 
Further, we assume fiL ^ oo, ^_r— >— cx). These assump- 
tions are necessary for the derivation of the Markovian 
dynamics of the array. 

Finally, we introduce a generic heat bath consisting of 
an infinite set of harmonic oscillators linearly coupled to 
the position of the central dot (Caldeira-Leggett model^^«) 
which simulates the dissipative interaction of the center- 
of-mass motion of the central dot with its environment 



Hhath + J^osc-bath + HqT — ^ ~ 



2 

\2mj 



W . 2-2 
W X 



(Id) 



The bath is characterized by its spectral density J(w) 



2 l^i 



■5[lo 



j). We take it in the Ohmic form^ 



2 ^-^j irij ujj ^ J } 

J{uj) = rn7a;/(^) where we have introduced the damp- 
ing coefficient 7 and /(^) is a model specific cut-off 
function /(x ^ 0) — *■ 1. As long as the cut-off fre- 
quency is much bigger than the frequency of the oscillator 
(lUc ^ luq) f would only contribute to the renormaliza- 



j?, + Auj'^ with Aw^ 



tion of ujq UJQ 

_1 r°da;iM ^ _27 r--dufi^). Here, we have ex- 

plicitly included the standard counter-term Hct can- 
celling this renormalization so that the bath solely in- 
duces dissipation and the cut-off function can be replaced 
by unity. 
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B. Generalized Master Equation 

For the description of the model we use the language 
of quantum dissipative systems^i As the "system" (or 
"device") we take the electronic states of the dots in 
the array (including the unoccupied state) plus the one- 
dimensional oscillator describing the center-of-mass mo- 
tion of the central dot. The electronic leads coupled to 
the outer dots and the heat bath interacting with the 
center-of-mass degree of freedom of the central dot con- 
stitute the reservoirs. The Hamiltonian of the system is 
then Ho = Hose +^^ci + i?ci-osc- For further reference we 
also introduce the Hamiltonian of all mechanical degrees 
of freedom, i.e. of the oscillator and the bath, reading 
■^osc = ^osc + ^osc-bath + ^bath + ^CT- The task is now 
to integrate out the degrees of freedom of the reservoirs 
to end up with an equation of motion for the system den- 
sity operator. We outline how the derivation proceeds in 
two steps first integrating out the leads in the high bias 
limit and then the heat bath in the weak coupling limit to 
get a generalized master equation (GME) for the system 
density operator. 

As in previous papers;iSi22ii£ we work in the high bias 
limit in which the bias between the leads is much higher 
than any other involved energy scale but the charging 
energy (cf. Ref. and Fig. ^ . The high bias assump- 
tion together with the wide-band limit means that after 
integrating out the leads the resulting dynamics of the 
system and heat bath is still Markovian. Following the 
derivation by Gurvitz and Prager'*? one can obtain the 
equations of motion for the density matrices (t'^"^ {t) of 
the system plus heat bath resolved with respect to the 
number of electrons n which have tunneled to the right 
lead by time t. We use the block notation analogous to 
the one used in Ref. {h = 1 throughout the paper 
except for figures): 

^^0^ = -^[HL, & r^o^ + ra^nn'^ - - o, l, . . . 

a'j'}^ = -^{I\[Ho, + H'o^, + i/ei-osc, J) 
+ (/|/Cdriv<T("V> for I,J^L,C,R . 

(2) 

Here ajj = (/|(t|J) are still operators in the oscillator 
and bath space. The "driving" kernel /Cdriv due to the 
coupling to the leads acts non-trivially only on the elec- 
tronic degrees of freedom and as unity on all the others. 
Hence also it can be written in the block notation 

-W2 (3) 

-o'7?l/2 -(7_r,c/2 -a-jiji J 

where the tunneling density of states F describes the 
injection rate from/to the leads. We still have to con- 
sider the off-diagonal block elements of the density ma- 
trix aoi, a 10 with / = L,C\R. They describe coher- 
ences between system states containing different number 
of electrons. In the formalism by Gurvitz and Prager*^^ 



these off-diagonal elements are identically zero by the 
construction of the theory (see also Ref. ilGj) . In other 
works, e.g. in Ref. l46l they can in principle appear, at 
least indirectly. In any case, whatever method is applied 
to our system, they are always decoupled from the rest 
of the elements. Moreover, they do not enter any expres- 
sions for quantities of physical interest that we consider, 
and can therefore be discarded. 

The GME for a{t) = J^n ^''"^ is found by summing © 
over n with the boundary condition^^a ct^"^) = 0. Due to 
this boundary condition the GME for a{t) is formally the 
same as Q just with the superscript index (n) omitted. 
This GME is used in subscction llV Cl and Appendix IBl in 
the sequential tunneling limit to derive a rate equation, 
from which both current and noise can be calculated, and 
compared to the full numerical evaluation. 

In general, there is no simple approximative analytic 
treatment of the problem nor is a direct numerical solu- 
tion possible due to the presence of the infinite number of 
bath degrees of freedom as a part of the system. To pro- 
ceed we have to integrate out the bath degrees of freedom 
to be left with the electronic and oscillator degrees of free- 
dom only which can be handled numerically. This could 
in principle be done in the weak coupling limit between 
the device and the heat bath by a perturbation expan- 
sion in the Cj's. This would amount to finding the "free" 
evolution of the device first, i.e. the evolution without 
the coupling to the heat bath but with coupling to the 
leads included. However, this free evolution is not unitary 
which significantly hinders any attempt to proceed. Even 
in the case of small coupling F to the leads, when the driv- 
ing Liouvillean is neglectedSS, one should diagonalize the 
device Hamiltonian (including the electromechanical cou- 
pling) and use the exact eigenenergies and eigenvectors 
as the input into the weak coupling prescription's*^ as 
was recently done in a dissipative double-dot system in 
Ref.lU 

Rather than following this lengthy procedure, we 
used the standard quantum optical damping kernel 
for a single harmonic oscillator in the rotating wave 
approximation^^iS also used in previous studies iiS42Sii& 
Strictly speaking, this can be justified only in the case 
of weak electromechanical coupling and small injec- 
tion. Nevertheless, we believe that the genuine non- 
equilibrium phenomena described later on are captured 
qualitatively correctly even with this kernel since the ker- 
nel mostly serves just as a "convergence factor" to sta- 
bilize the stationary solution. As will be seen below, 
the sequential tunneling limit is extremely well captured 
within the adopted approach. This is perhaps not too 
surprising since in that limit the coherence between dif- 
ferent dots is negligible. On the other hand, the clear 
advantage of our choice of the damping kernel is that it 
preserves charge conservation throughout the whole cir- 
cuit while this may not happen in general in the weak 
coupling prescription (see Section [HI D|) . Refinements of 
the present approaches to deal with the above issues are 
in our opinion a challenging task for the future modelling 
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of NEMS. We would like to point out that the above men- 
tioned concerns about additivity of the two baths apply 
also to the case of the one-dot setup traditionally used 
for the description of the shuttling phenoniena22iSi2i but 
the problem stemming from the coherence present within 
the array is absent there. 

Bearing all these cautions in mind, we are ready to 
state the generalized master equationi^ for the n-resolved 
density matrix of the system: 

POO — ~*[-nosc, Poo J + -i-dampPoO ~ PoO ~^ ^ PrR 
'PU = -*(/|[i?cl + i?osc + i?el-osc,p("^]| J) 

+ AampP^'}' + (/|/:drivp("' I J) for I,J=L,C,R. 

(4) 

The commutator terms in (@J) describe the coherent evo- 
lution of the isolated device. The driving kernel £driv is 
given just by substitution a ^ p in J^J: 

(Poo -pLfl/2\ 

-pcr/2 . (5) 

-Prl/2 ~prc/2 ~prr J 

Finally, the damping kernel^^ (acting as unity on the 
electronic degrees of freedom) reads 

>CdampP = -■^n(aa'^p - 2a^ pa + paa'') 

, (6) 

— — (n + l){d' ap — 2ap<r + pa' a) 

where 7 is the damping rate and h — ns(wo) = 
(exp(cL;o/^_BT) — is the mean occupation number of 
the oscillator at temperature T. This term describes the 
effect of the environment on the oscillator, consisting in 
mechanical damping and random quantum and thermal 
excitation (Langevin force). The issue of the appropri- 
ate choice of the damping kernel is, however, quite subtle 
in many respects even in the case of a simple harmonic 
oscillator used here. There is a well-known dilemma be- 
tween the rotating wave approximation form (conserving 
the positive definiteness of the resulting density matrix) 
which we use in this work versus the translationally in- 
variant form (yielding correct equations of motion for the 
mean coordinate and momentum) used previously.^^'"^^- It 
is known that this dilemma cannot be solved within the 
Markov approximation (without relaxing the condition 
of approach to the canonical thermal equilibrium state 
for asymptotic times; for a thorough discussion of this 
issue see Ref . 153^) . We have carried out a number of nu- 
merical checks, and have found out that in the present 
case there are only minor differences in the obtained re- 
sults. A practical advantage of the present choice is that 
it leads to faster numerical convergence. 

We can recast the GME into a compact form 

/5(") = (£-Xofl)p<"^+W""'' , 

p = £p with P = Y, P^"^ aiid =0 

n=0 



where Tqrp = r|0)(i?|p|i?)(0| (the symbol Tqr denotes 
the superoperator of the particle current across the junc- 
tion OR between the right dot and the right lead, for a 
discussion on superoperators see below). 

The dynamics of the device described by the above 
generalized master equation (UJ constitutes a quantum 
Markov processm^ The Liouvillean C determines the evo- 
lution superoperator exp{Ct) which fully characterizes 
the resulting quantum Markov process. It can be used 
to calculate arbitrary multi-time correlation functions of 
any system operators, i.e. operators acting as unity on the 
Hilbert space of the reservoirs, by using the multi-time 
structure of the quantum Markov process (often referred 
to as the quantum regression theorem) — for details see 
Ref.El Sec. 5.2 or Ref.ls^, Sec. 3.2. Therefore, not only 
the mean value of the stationary current within the array 
as in Refs. 0, 113 can be evaluated in this way, but also 
its higher order correlation functions, in particular the 
current noise spectrum, become accessible. The calcula- 
tion can only be done for the junctions within the array. 
For the outer junctions between the outer dots and leads 
the quantum regression theorem cannot be applied since 
the corresponding current operators involve the lead elec- 
trons, thereby not being system operators. However, the 
n-resolved form of the GME iQ enables us to calculate 
the current noise spectrum also for those junctions. Both 
methods yield equivalent results as we will show later in 
Section UTTDl 



C. Notational details 

The linear operator C which acts on the density op- 
erators, as specified by (©-0, can be handled (at least 
formally) as any other linear operator. We can asso- 
ciate a matrix (infinite in our case) with it and perform 
standard linear algebra operations. In order to avoid 
confusion with "normal" quantum mechanical operators 
acting in the "normal" Hilbert space of the system, the 
vector space of "normal" operators is called the Liou- 
ville space or the superspace, and the Liouvillean and 
other linear operators acting in the superspace are called 
superoperators (or supermatrices). In the following, all 
superoperators will be denoted by calligraphic symbols 
and the vectors of the superspace in the bra-ket nota- 
tion will be distinguished from the normal vectors in the 
Hilbert space by double brackets, e.g. V ^ \v)) with V 
being a "normal" quantum mechanical operator. 

If is an orthonormal basis in the Hilbert 

space of the system then all the projectors {|r7i)(n| = 
l'^'^))}mn=i form an orthonormal basis of the corre- 
sponding Liouville space with respect to the scalar prod- 
uct ((a 1 6)) = Tisysi^^B). The matrix representation of 
superoperators follows analogously to the normal Hilbert 
space case, i.e. O = I]fei,Tnn I^O)((^^|C'|'Tm))((TOn| = 
Sfei.mn 1^0) Oki,mn{{mn\. There is a unique mapping be- 
tween matrices representing the operators in the Hilbert 
space and the vectors in the Liouville space. The 
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operator O = J2k i\^)Oki{l\ represented by the ma- 

/Oll Ol2 ... \ 

O21 O22 ... 



trix 



corresponds to the vector |o)) = 



J2ki^ki\kl)) represented by the column vector o — 
{On, Oi2,Oi3, ... ,021,022,023, ■■ - V (the exact order- 
ing depends on the chosen ordering of the double indices 
kl). Therefore, we will in the following use the two rep- 
resentations interchangeably. 



III. NOISE CALCULATION 

A. Definition and properties of the current noise 
spectrum 

In this subsection we define the current noise spectra 
for different junctions present in our model and analyze 
several of their properties. First, we find the current op- 
erators across different junctions. From the equations of 
motion for the operators of the occupation of the respec- 
tive dots hj = \J){J\, J — 0, L,C, R reading 

e-jj nj = -ie[nj, H] = /,/+ - /j_ (8) 
at 

we identify the corresponding charge current operators 
(electronic charge is e < 0; electrons flow from left to 
right) 



dt 



CkL 





k 




h- 


= ic+ = 


icL = 


ic- 




Irc = 


Ir- 




Iqr = e 



d_ 

di 



le 



j2VkR{\R){o\ckR-ci^mR\ 



(9a) 

(9b) 
(9c) 

(9d) 



with Nl = J2k^kL^kL, Nr = J^k^kR^kR being the op- 
erators of the number of particles in the left and right 
lead, respectively. 

We next define different current-current correlation 
functions (a, b = LO, CL, EC, OR) 

Cab{T)= lim \h{i^{t + T),h{t)}) -{ia{t + T)){ii,{t)) 

t^oo L Z 

= hm i({A4(i + T),A/fc(t)}) , 

t—*oo Z 

with Aiait) = iait) - {ia{t)) 

(10) 

which in the stationary limit are functions of r only. We 
also note the property Cab(— t) — Cbaij). The current 
noise spectrum ism 



(11) 



The diagonal elements Saai^^) of the noise matrix are 
non-negative as can be shown by using the Lehmann rep- 
resentation. 

In general, for an arbitrary frequency the noise depends 
on the position where the current is measured. How- 
ever, in the limit cj — > charge conservation implies that 
the noise becomes independent of the measurement po- 
sition along the circuit, i.e. SaaiO) = Sbb{0) = S'at(O) — 
Sba{Q), a^h and it also equals the shot noise component 
of the spectrum measured in the leads. This statement 
is proven by considering current correlation functions for 
two adjacent junctions J-|-, J—S^ The charge conserva- 
tion condition JSjl gives 



C. 



j+j- 



_(r) = i({A/,+ (r), A/^+}) = i({A/V(r), A/,,+ }) 
+ ~({eAhjir),AIj+}) 



= Cj^Mr) + i^({eAn,(r), A/^+}) 



(12) 



which implies S'.7+j+(0) = 5'j_j+(0). The relation 
Cah(-T) = Cbair) yields Sj-j+{-uj) = Sj+j-{uj) and 
by using the charge conservation again we can finally es- 
tablish Sj+j^{0) = Sj-j-{0). Altogether we find that 
the zero-frequency noise is the same for any combina- 
tion of the junctions, i.e. S'ab(O) — S{0) > for any a, b 
(not necessarily adjacent; this generalization is straight- 
forward). 

The current operators Ic l , Irc H9b|l , (Pc|l between the 
dots are obviously system operators in the sense that they 
operate as unity on the degrees of freedom of the leads 
and the heat bath. Therefore, we can use the formal- 
ism of quantum Markov processes to evaluate correlation 
functions involving these operators using the quantum re- 
gression theorem — this will be done in subsection IIIIBI 
This is not the case for the operators of current between 
the outer dots and leads Ilo, Iqr given by (1^ . H9d|l . 
However, the noise spectra across these two junctions can 
still be calculated using the n-resolved form of the GME 
with the help of the following identity for the zero- 
frequency current noise (for the junction OR, the case LO 
is analogous) 



d 
di 



{{Qlit)) - (Qfl(t)>- 



dTCQR,QR(T) = SoR.ORiO) 



with Qnit) = eNBit) - eiVfl(O) = / dt'ioR{t') 

Jo 



(13) 



This identity suggests the interpretation of the zero- 
frequency current noise as the "charge diffusion 
coefficient"— and will be used in subsection IIII CI for an 
alternative evaluation of the zero- frequency current noise. 
The equivalence of the two approaches is shown explicitly 
in subsection IIII DI 

We finally comment on the physical relevance of the 
noise spectra calculated in this paper. Since the zero- 
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frequency noise is position-independent the noise calcu- 
lated for the junctions within the system should also be 
measured in the leads. However, in practice there is al- 
ways the important 1// contribution to the noise which 
actually dominates experiments for very low frequencies 
and which is not accounted for in our model. Therefore, 
as mentioned in Ref. the measurements of the shot 
noise must be performed at non-zero frequencies of the 
order of 1 kHz where the 1// noise component becomes 
insignificant. However, the shot noise measured in this 
way is still appropriately described by the zero- frequency 
current noise calculations since its typical frequency scale 
is of the order of 1 THz. 



1 — V. Using these projectors and the relations VC = 
CV = 0, £ = QCQ the resolvent can be expressed as 

g{-iuj) = i-iuj - C)-^ = {-iLuV - iujQ - QCQ)'^ 

lUJ ILO -\- L 

K. — ^-V - QC^^Q = — ^-V - n for small oj 

iuj iuj 

(17) 

where we have defined the pseudoinverse of the Liouvil- 
lean TZ = QCt^Q. Substituting the term — ^ in the first 
term of H16|l gives 



B. Quantum regression theorem (QRT) 

With QRT it is possible to calculate the current noise 
within the system (i.e. for Icl, Irc)- For r > QRT 
gives (cf. Ref.m Sec. 5.2) 



1 

2iLU 



Tr3ys(/a|0))((0|{/b,p^*^n) 



Cabir) = 77Tr,y,(4 eM^r){h, p''*^*}) - (14) 



for a,b = CL,RC, where / = limt^oo(-^a(i)) — 
Trsys(/aP'^'''*) is the stationary current (constant through- 
out the circuit). In case r < we use the symmetry 
property Ca&(— r) = Chair). Now, let us evaluate the 
spectrum 



Sab{i^) 



dTCab{T)e'' 



dTCab{T)e^'^^ + / dTCba{r)e 
Jo 



(15) 



We consider in detail the first term denoted S'^{,(ci;), the 
second one {Sf^^{Lo)) follows analogously. Introducing a 
convergence factor cj — > cj + iO we get 

S+,{oj) = iTr,y,(4(-^c.-£)-i{4,P^*^'}) + ^^'- (16) 

Since we are interested in the limit cj — > in the end 
we have to handle somehow the singularities associated 
with the resolvent Qi^—ioj) — {—iuj — C)^^ and the sec- 
ond term in Hlti|) in that limit. The problem with the 
inverse of C is the existence of the unique null vector 
|0)) which is proportional to the stationary density ma- 
trix because Cp^'^^* = 0. There exists a corresponding left 
eigenvector belonging to the zero eigenvalue of C denoted 
by ((0| which is not just the hermitian conjugate of |0)) 
(i.e. ((0| 7^ |0))^) because £ is non-hermitian. However, 
since Tr, 

that ((0| ^ 1, i.e. {{0\C\a)) = Tr 

Thus, we have |0)) ^ p"*'^*, {(0\ ^ 1 with ((0|0))^= 1 al- 
lowing us to define the projectors V = V^ = |0))((0|, Q - 



sys{CA) — for any system operator A we deduce 

.sys{iCA)=0. 



--^Tr,y,(/,p^*^')Tr,y,({/,,p^''^'}) = 



(18) 



which cancels the last term of (|16|) . Applying the same 
procedure to 5*^(0) we find 

5,,(0) = 5+(0) + ^-(0) 

= -i(Tr,y,(47^{/6, ff''^'}) + Tr,y,(/^,7^{4, ff^'^^})). 

(19) 

If we introduce the superoperators of (particle) current 
TcL, Irc defined by their action on the system density 
matrix as follows eXaP — ^{Imp}, a = CL,RC with the 
property / = eTisys^aP^^'^*' = e((0|la|0)) we can rewrite 
the above equation in a compact form 

Sab{0) = -e'miaKIb + Xfc7^X,|0)) a, b = CL, RC . 

(20) 

This equation constitutes the main formal result of this 
subsection and forms the basis for further formal manip- 
ulations and eventually the numerical treatment. 



C. Counting variable approach — evaluation of the 
charge diffusion coefficient 

Using the n-resolved form of the GME lO we could 
in principle find the full counting statistics (FCS) of the 
charge transfer through the junction between the right 
dot and the right lead, i.e. the probabilities P„(i) that n 
electrons tunneled into the right lead across the junction 
by time t given by Pn{t) = Trsys/5^"-'(<). Here, we are only 
interested in the evaluation of the zero-frequency noise 
for which we just need the mean and the mean square 
charge tunneled into the right lead by time t given by 

(QflW) = eE„"^n(*), {QUt)) = e^En^^Pnit). Us- 

ing the definition of the current IjQdp and the identity 
(|13|l we find the stationary mean current and the zero- 
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frequency current noisei^ 

IoR = e^y2nPnit)\ =eVnP„(t)| , (21) 
at — U— *oo ^ — ^ U— foo 

n n 

d ' 



dt 



Because the \aige-t behavior of F{t; z) is related to the 

small-s behavior of F{s\ z) we study the asymptotics of 
the above expressions as s — > 0+. This is entirely deter- 
mined by the resolvent Q{s) in the small-s limit. We can 
use the results from the previous subsection and substi- 
tute —iLO s to get the leading asymptotics of G{s) for 
s 0+. Thus, we obtain 



t — >oo 

(22) 

We evaluate Pn{t) from Eq. and find 

Pnit) = Tr,ys[loR{p^''-'\t) - p(")(t))] (23) 
and consequently 

n 

Y^^Pnit) = Tr,y,(XoflE/5(")(t)) = Tr,y,(Xofl/5(t)) , 

n n 

(25) 



1 



d ~ 
az 



(32) 



7ei-oi?7'p(0)-7'En/5(")(0)j . (33) 



In the time domain this gives 
F(t;l)|*^oo 



d_ 

dz 



F{t-z) 



z—l,t—*oo 



^stat 



(j 



t + C 



init 



(34) 
(35) 



Yii?Pn{t) = Tr.ys Xoi?,(2E?7p(")(i) + p(t)) , (26) where C™* = Tr,y,( ^„ 7ip(") (0) - Xoi^7^/5(0)) is an ini- 



where according to the definition X]n/'^"''(0 = P(0- 
Now, we employ an operator-valued generalization of 
the standard generating function technique to calcu- 
late X]ri"/''"H*)- We introduce the object F{t;z) — 
^^/5(")(t)z" which has the properties F{t;l) ~ 

p{t), -g^F{t; z)\z=i = Sn'^p'"^*) ^^^d satisfies the equa- 
tion of motion 

^^F{t-z) = {C+{z-l)IoR)F{t;z) . (27) 

Using the generating function the current noise formula 
can be rewritten as 



Fit- 1) 



tial conditions dependent constant and the stationary 
current is given by / = eTrsys(Ioi?/5*'*^')- The corrections 
to the large time asymptotic behavior are exponentially 
small — the approach to the stationary state in a Marko- 
vian system is exponential. In particular, it is important 
that there is no j correction to F{t] l)|t^oo (which would 
correspond to a In s-like divergence in the resolvent as 
s 0+) since it would combine with the linearly in t di- 
a 



vergent term in -0- F(t; z) 



Sor,or{0) = (^Trsys Ioii(2-^F{t; z] 
- 2Tr,y, (loRpit; 1)) Tr,y, (g^Fit; z) 



to yield a finite term 

z — l^t — »oo 

in if^ . We substitute the above asymptotic formulas 
into Eq. (|28|l . use the definition of the stationary current 
and the identities Trsys/5''*''* = 1, Trsys7?.» = to get the 
final result for the zero- frequency current noise at the OR 
junction. 



(0) 



el - 2e^TTsys{loRT^1oRP'' 
e2((0|Xofl-2Joi^7^Xoi^|0)) 



(36) 



(28) 



The equation of motion for F{t; z) H27|l can be solved via 
the Laplace transform F{s; z) = dte~''^F{t; z) giving 

(s-/:-(z-l)Xofl)F(s;z)=Ep(")(0)z" , (29) 

n 

with p^")(0) being the initial conditions. Recalling the 
definition of the resolvent Q{s) = [s — C)^^ of the Liou- 
villean we arrive at 



F{s-A)^g{s)m 



(30) 



— ~F{s- z) = g{s)IoRg{s)m + e(s) V (0) . 

OZ z^l ^ — ' 

n 

(31) 



In the algebra leading to H36II the linearly divergent 
terms in t and the initial condition terms cancel identi- 
cally so that we are left with a regular, initial-condition- 
independent expression as expected and necessary. Sim- 
ilarly, for the LO junction one finds 

^LO,Lo(0) = e/-2e2Tr,y,(XLo7^XLO/5^'^*) , , 
' ~ 37) 
= e2((0|XLO"2Xio7^XLo|0)) 

withXio/5 = r|L)(0|/)|0)(L|. 



D. Equivalence of different approaches 

We show the equality between the expressions (|2(JII and 
(|36|l . (|37|l . Both formulas contain the same basic building 
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block consisting of terms of the type laTiTb- However, 
there is an obvious difference: the presence of the so- 
called self-correlation or Schottky term (proportional to 
the mean current) in formulas Ij^GI) . (|37|) . Yet, they give 
the same value for the zero-frequency noise in the end as 
we now proceed to show. 

The independence of the zero-frequency noise from the 
position along the circuit has been shown quite generally 
in subsection IIII Al using the charge conservation. Thus, 
the only task now is to find the corresponding expression 
for the charge conservation within the superoperator lan- 
guage. Following the purely stochastic analogy^® we find 
that the charge conservation condition JSJ is expressed 
in terms of superoperators by the following equation 

[Nj,C]=Ij+-Ij- (38) 

with the superoperators of occupation of the "site" 
J, J = 0, L, C, R being given by Afjp — ^{\J){J\, p}, the 
current superoperators were defined previously and 
the convention for J± is the same as in Eqs. 0. The 
above relation follows from the definitions of the respec- 
tive quantities and equations Q)-®. 

Since the heat bath does not couple directly to the 
electronic degrees of freedom its degrees of freedom do 
not enter explicitly the current and occupation operators, 
cf. ijSJ and and are therefore absent from the corre- 
sponding superoperators. We believe that this property 
should be reflected in the identity [N'j,Cdamp] = for 
any choice of the damping kernel. Obviously, this condi- 
tion is fulfilled for our choice of the damping kernel 
However, for the generic weak coupling prescription'*^*^ 
for the damping kernel the above identity may not be 
satisfied which would break the charge conservationiS 
This raises the possibility that there is another problem 
with the Markovian weak damping prescription analo- 
gous to the translational invariance issue threatening the 
charge conservation for damped NEMS involving coher- 
ent charge transfer (such as our quantum dot array). 
This issue deserves further investigation. 

The charge conservation relation H38|) is used to prove 
the position independence of the mean current / ~ 
e((0|Xa|0)) and the zero-frequency noise Sab{0) for any 
a, b. The mean current conservation follows from 

/ = e((0|X;+|0)) =e((0|X,_|0))+e((0|[A(7,£]|0)) 
= e((0|Xj_|0)) due to C\0)) = 0, ((0|£ = . 

Analogously, we prove the equivalence, for example, be- 
tween S'oj?,oi?(0) (EH) and Srcm^c{0) llOll- Substituting 
for J = i? into the expression (pn|l for S'i?c,flc(0) we 
get in the first step 

SRCRciO) = -2e^{b\lRcniRcm 

= -e\{6\ioRniRc +iRcmQR\o)) 

= SRcoRiO) — SoR,RciO) 

(40) 



bearing in mind CTl = TZC = Q = 1 — |0)) ((0| and finding 

e[XRc,AfR]p = liilRC, \R){R\],p] which yields zero when 
traced over. We proceed similarly in the second step and 
obtain 

SoBmiO) = -2e2((6|Xoi^7^Xoi^J0)) +e2((0|[Xofi,AAfl,]|0)) . 

(41) 

The second term can be evaluated as [Xqr^JVr] = 
[JVojTqr] = Tqr recovering finally the expression (|36|l for 

'S'o_R,Ofl(0). 

By extending the argument to other combinations of 
the junctions we can summarize the formulas for the 
zero-frequency noise 5(0) — S'/+_j+(0) for any I,J — 
0, L,C, Rin the compact form as (compare with the anal- 
ogous expression for the purely stochastic case in Ref.l56l 
Eq. (26)) 

s{0) = -e'mii+nij+ + i.,+nii+ |o)) 

+ W((0|[A/},Xj+]|0)) for any/, J. 

This equation merges the two approaches into a single 
picture unifying both the pure quantum mechanical and 
pure classical stochastic formalisms. It has a quantum- 
mechanical-like form of a "mean value" of the pseudoin- 
verse of the Liouvillean symmetrically flanked by two 
current superoperators corrected with the classical-like 
self-correlation term. The self-correlation term is only 
effective for the diagonal elements of the current-current 
correlation matrix and, moreover, is non-zero just for the 
outer junctions where it contributes by the mean current. 

E. Notes on numerical evaluation 

From the results obtained thus far we see that the eval- 
uation of the noise involves two steps. At the first step 
we find the stationary state ff'^'^^ = lim(_»oo exp(£t)po 
independent of the initial condition po and equivalently 
given by the equation 

^^stat ^ Q ^ Trsysp'*^' = 1 . (43) 

Having found p***^* we can fully characterize all one-time 
quantities pertaining to the system such as occupations 
of the different dots, mean current, Wigner functions of 
the oscillator in different charge states etc. 

To evaluate the noise (second step) we have to find 
the pseudoinverse of the Liouvillean TZ — QC^^Q. In 
practice, we actually do not have to evaluate the whole 
pseudoinverse but we fix a given combination of junctions 
and evaluate the auxiliary quantities Ea = eRXaP^^'^^ de- 
termined by the equation 

eta = eX.p^'-* - Iff'-' , Tr,y,Ea = . (44) 

Eq. H44|l has a solution since the right hand side lies in 
the range of C (the trace of the right hand side is zero) 
and the freedom of adding any multiple of the null vec- 
tor to a particular solution is fixed uniquely by the trace 
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condition TrgygSa = 0. Of course, this is equivalent to 
the uniqueness and regularity of the standalone pseudoin- 
verse TZ. Moreover, TZ preserves hermiticity so that the 
quantities Sa are Hcrmitian as they should be to give a 
real zero-frequency noise. This follows from the property 
(CA)'' = CA for any Hermitian A and the trace- fixing 
condition Tv^yJ^a = of Eq. ll^ . 

Equations ()43|l and H44() form the starting point for 
the numerical implementation of the noise calculation. 
After the truncation of the oscillator Hilbert space to 
the N lowest energy states the size of the supermatrix 
C becomesS^ lOiV^ x lOA''^ which makes direct calcula- 
tions prohibitive due to memory and computation time 
requirements for any realistic N of the order of 30-40. 
These problems with the numerical implementation of 
the superoperator techniques can be circumvented by 
employing iterative methods in which only the proce- 
dure/routine yielding CA for a given A is needed.-—' Ob- 
viously, this does not require the storage of the whole 
supermatrix C. On the other hand, as with any itera- 
tive method, the convergence of the iteration becomes 
an issue. In Appendix ^ we give a brief review of the 
usage of the Arnold! iteration in our calculations. Its in- 
tent is to guide the reader through the algorithm so that 
it can be reproduced with the help of the mathematical 
references i^SiS 



IV. RESULTS 

We now turn to the numerical results for the mean cur- 
rent /, zero-frequency noise S'(O) — Sah{Q) (for any a,b 
— see above) , and the Fano factor F — as functions 
of the device bias £(, for different sets of the other param- 
eters. First we present a generic plot in the parameter 
regime considered by Armour and MacKinnon and com- 
ment on the general features which we can observe in it. 
We then give a tentative interpretation of those features 
supported by phenomenological arguments and results 
found in different limiting cases studied further on. In 
particular, we consider two specific limiting cases where 
at least a partial comparison with approximate analytic 
theories can be made, namely (i) the limit of small damp- 
ing which is relevant for the issue of shuttling and strong 
inelastic cotunneling and (ii) the limit of weak inter-dot 
coupling which implies in a certain device bias range the 
sequential tunneling regime. 



A. Generic case 

In Fig. 121 we plot the mean current, zero- frequency 
noise and the Fano factor as functions of the device bias 
and temperature for one of the parameter sets considered 
in Ref. 16. We include non-zero temperature and extend 
the device bias range considered previouslyi^ to negative 
values which is relevant for non-zero temperature. 



The dotted lines show the results for the static ar- 
ray. By applying the theory of Sec. Illll to the static 
array we found analytic expressions for both the mean 
current and the Fano factor which wc, however, do not 
present explicitly here since the formulas are quite in- 
volved. The mean current has a resonant peakSS around 
Eh = while there is a dip in the noise around Sb = 
which was also found analytically for a two dot array by 
Elattari and Gurvitzi^ They attributed the dip to the 
strong Coulomb interaction on the array. Our Fano fac- 
tor shows a crossover from the sub-Poissonian (F < 1) 
dip around £fc = to super-Poissonian (F > 1) "shoul- 
ders" starting around £b ~ zbVbe""^" which approach 
the Poissonian limit F = 1 for large device bias. The 
Poissonian limit of the Fano factor for large Sb is un- 
derstood when one notices that the current in that limit 
is very small. Therefore, electrons tunnel through the 
array sparsely and, consequently, there is no correlation 
between successive tunneling events which form a clas- 
sical Poisson process with the (Poissonian) value of the 
Fano factor F = 1. While the dip around zero and the 
Poissonian limit for large device bias were observed in 
the two dot case as well^ the Fano factor exceeding one 
was not present there. We attribute the super-Poissonian 
behavior to the (elastic) cotunneling through the central 
dot. 

Now, let us discuss the results for movable arrays. The 
characteristic features are the peaks in current and noise 
at the device bias around a non-zero integer multiple of 
the oscillator frequency due to electromechanical reso- 
nances. The current peaks at zero temperature (there- 
fore, only for positive multiples of the frequency) were al- 
ready observed in previous works. ^^'^^ Some of the noise 
peaks have further fine structure which is even amplified 
in the Fano factor exhibiting a rather complex behav- 
ior around the peaks, especially at low temperature, and 
showing also strong temperature dependence. 

The zero device bias behavior is clearly governed by 
the static array physics which is due to partial decou- 
pling of the electronic and oscillator degrees of freedom 
at £{, = when the electrostatic interaction on the cen- 
tral site —^^x\C){C\ is turned off. The remaining inter- 
action stemming from the x-dependence of the hopping 
amplitudes ti^{x), tji{x) is too weak to modify the static 
result in the vicinity of eb = even for high temper- 
atures. Some discrepancy between the static and high 
temperature dynamic cases around = is found for 
higher values of a ~ 1 (strictly quantum case from the 
oscillator point of view which was previously studied in 
the one-dot shuttling sctup^^ '^''), yet the effect is not very 
pronounced anyway (not shown). 

The peaks at non-zero multiples of the oscillator fre- 
quency were already previously attributed to electrome- 
chanical resonances liSiSS Yet, this explanation is rather 
broad and covers a range of processes which can be re- 
sponsible for the electronic transport such as cotunnel- 
ing, phonon-assisted tunneling, or shuttling occurring 
around different resonance peaks^^i^ The discrimina- 
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FIG. 2: The mean current I, zero-frequency noise ^(O), and the Fano factor F as functions of the device bias et for the 
static three dot array (dotted line) and the vibrating array at different temperatures given by the mean oscillator occupation 
number n. The other parameters are Vo = 0.5hLoo, a = 0.2V2^mujo/ti, xo = {5/V2)^h/mwo, 7 = 0.025ljo, T = O.OStJo which 
corresponds to the case studied in Ref. Ilal . Fig. 6. 



tion between the dilFerent processes is quite complicated 
since it cannot be inferred directly from a single / vs. Sb 
curve. Either one has to study the dependence of the 
curves on different parametersi^ or some other kind of 
information about the system must be obtained. A pow- 
erful choice is to calculate and analyze the Wigner dis- 
tribution functions of the oscillator in the phase space 
(possibly charge-resolved) These characterize the 

state of the system very well and we will use them in this 
study too. However, even though they are an excellent 
theoretical tool to study NEMS their connection to data 
extractable from a real NEMS experiment is at best re- 
mote. Therefore, diagnostics based on the measurement 
of the current statistics is clearly preferable and, there- 
fore, our aim is to correlate particular features observed 
in the noise with specific transport mechanisms within 
the array as identified by the theoretical analysis involv- 
ing also phase space plots. 

To achieve this goal we will study different limiting 
cases in which particular features of the noise (more pre- 
cisely of the Fano factor) are pronounced so that they can 
be attributed to specific transport mechanisms. Yet, the 
results do not allow to associate a given value of the Fano 
factor to a specific mechanism. It is more reading of the 
whole / vs. Sb curve at least locally around a peak which 
gives us the notion of what mechanism(s) are involved in 



the transport at that given peak. 

As a rule of thumb we can say that the super- 
Poissonian peaks of the Fano factor correspond to cotun- 
neling through the central dot. This statement is sup- 
ported by the limiting studies discussed below, and also 
by the following evidence from Fig.|21 The peaks only oc- 
cur for small temperature and disappear with its increase 
pointing out to a coherent effect. They also appear pre- 
dominantly at odd multiples of the oscillator frequency 
which is consistent with the cotunneling picture between 
the outer dots excluding the central one due to the en- 
ergy mismatch. On the other hand, the dips in the Fano 
factor curves are due to some form of the sequential tun- 
neling via the central dot. The most important aspect is 
that the process proceeds via a real intermediate state on 
the central dot in contrast to the virtual nature of the co- 
tunneling process. The real sequential process is subject 
to the charge conservation which is a strict law strongly 
suppressing the Fano factor— and causing the dip. The 
sequential tunneling picture still involves different mech- 
anisms distinguished by the detailed state of the oscilla- 
tor. The oscillator might be in a general non-equilibrium 
state during the tunneling events (this scenario encom- 
passes both the shuttling^ and a general non-equilibrium 
oscillator-assisted tunnelingS^ mechanisms) or it could 
equilibrate between consecutive tunneling events. The 
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latter case is studied in detail in subsection IIV CI 

The two charge transfer mechanisms (cotunneling and 
sequential tunneling) may coexist, i.e. part of the cur- 
rent is carried by the cotunneling mechanism and the 
other part by the sequential tunneling, and their relative 
weights depend strongly on the parameters. For exam- 
ple, the transport around « 2fkvo is typically governed 
by shuttling which results in the dip while cotunneling is 
dominant around Sb ~ Shujo giving a peak. However, the 
dip around Sb ~ in Fig. |21 changes into a clear peak 
when a is enlarged up to a « 0.4 (not shown). This 
behavior is still not well understood. Even more com- 
plicated is the behavior around £{, « 'Ihioo where there 
is a dip in the peak. As we show in the next subsection 
this corresponds to a fast crossover between the cotun- 
neling and shuttling transport mechanisms in the vicinity 
of £{, = 4. In order to support the above statements for 
the generic parameters we study particular limiting cases 
which enable us to associate specific features of the Fano 
factor curves to specific mechanisms. 



B. Small damping: shuttling and strong inelastic 
cotunneling 

In this section results for small damping case, i.e. 7 ^ f 
with / a representative value of the current (given, e.g. by 
its value at the zero device bias peak), are presented^4 
First, we focus on the device bias range Sb ~ — 2.5fiwo 
where electromechanical instabilities which can be re- 
lated to shuttling were inferred indirectly from the be- 
havior of the mean current predicted by quasiclassi- 
cal studies, and subsequently directly observed in the 
phase space. The intuition and simple theoretical esti- 



mates (the zero-frequency noise is given by the ratio of 
the variance and the square mean of the waiting time be- 
tween consecutive loading events of the classical shuttle, 
see Eq. (4.48) in Ref. ^ suggest that shuttling is a low 
noise phenomenon with the Fano factor close to zero in 
the nearly perfectly developed shuttling regime. This was 
recently confirmed by more sophisticated calculations for 
the classical driven^^ and quantum'^^ shuttle in the one- 
dot setup. In the present, more complicated setup the 
shuttling is obscured by competing mechanisms (coher- 
ence between dots, strong Coulomb blockade affecting 
the whole array) and we will study the consequence of 
this fact on the behavior of the Fano factor. 

In Fig.|31we show the results for the mean current and 
the Fano factor for zero temperature and three different 
(small) values of the damping. In Ref. we presented 
the phase space plots of the oscillator which we introduce 
here in more detail later on (see Eq. I|45(l and Fig. inj. 
They described a similar parameter range and showed 
gradually developing shuttling around Sb ~ ^jJq, 2hwQ 
with increasing injection rate F. At these resonance 
points the current has peaks moderately changing with 
the increase of the damping and the Fano factor has lo- 
cal minima with possible shoulder-like structure further 
from the resonance points in case of the smallest damp- 
ing. As established more explicitly below, the shoulders 
are a signature of coherent processes through the whole 
array (cotunneling) and, therefore, are destroyed by the 
increased damping. 

At the same time the absolute values of the local 
minima of the Fano factor at the resonances become 
deeper by the increased damping. We conjecture that 
this somewhat surprising behavior can also be attributed 
to the destruction of the quantum coherence and to the 
crossover into the non-equilibrium sequential tunneling 
regime partially encompassing shuttling. The minimum 
of the Fano factor curve starts to increase again with fur- 
ther increase of damping (not shown) as expected from 
the classical shuttling theory. The minimal value of the 
Fano factor achieved for the given set of parameters was 
^rnin ~ 0.25 which corresponds to a partially developed 
shuttling regime and was also confirmed by the phase 
space pictures (not shown). 

Next, we focus on the range Sb ~ 2.5hujo — -i.bhLUQ in- 
volving two current peaks around Sb ~ 3fiwo, 4fia;o- As 
we already mentioned in the generic case the peak around 
£b « 3hoJo corresponds to cotunneling while the behavior 
around et « ifiiuo is given by a complicated interplay be- 
tween both mechanisms (cotunneling and sequential tun- 
neling). With lower damping the differences in the Fano 
factors of the two mechanisms become more pronounced 
as we show in Figs. 01 and [3 In Fig. 0] the mean cur- 
rent and the Fano factor as functions of the device bias 
£b are depicted for several (small) values of the damp- 
ing. We see the strong damping dependence of the mean 
current and the Fano factor around Sb ~ Shiuo and in 
the "shoulder region" around Sb « itujjQ. On the other 
hand the mean current as well as the Fano factor do not 
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FIG. 4: The mean current and Fano factor for Vb ~ 
0.76huJo,a = 0.28y^mujo/h,xo = 5y^h/mLuo,r = 0.2ujo,T = 
and different values of the damping coefficient (in units of 
ujo) in the strong inelastic cotunneling/shuttling regime. The 
dots on the curves corresponding to 7 = 0.0125 denote the 
points for which the Wigner functions in Fig.|Slare plotted. 




FIG. 5: Phase space representation of the oscillator around 
the transition from the shuttling to the strong inelastic co- 
tunneling regime at Eb/hiJo = 3.96,4.04,4.12, respectively 
(columns from the left to the right). The respective rows show 
the Wigner distribution functions for the empty (Wuu) or oc- 
cupied (Wcc) central dot, and the sum of the two (Wtot = 
Wuu + Wcc) in the oscillator phase space (horizontal axis - 
coordinate in units of Ayft/rntJo, vertical axis - momentum in 
\/hmujo, the grid is at 2.5 in the dimensionless units). The 
other parameters are: Vb = 0.76fttJo,a ~ 0.28y' muo/h, xq ~ 
5-^/l/mtJo, 7 = 0.0125tJo,r = 0.2uJo,T = 0. The parameters 
correspond to the dots in Fig. |1] The Wigner functions are 
normalized within each column. 



depend strongly on the damping in the close vicinity of 

We attribute the first type of behavior to cotunneling. 
It is manifested by a strong damping dependence of the 
current and the Fano factor, the Fano factor reaches very 
high values of the order of F « 50 for small enough damp- 
ing. The threshold for the quasi-divergent behavior of the 
Fano factor is roughly 7thrcsh ~ f ; for the damping befow 
this threshold the Fano factor starts to increase. We want 
to point out that a giant (divergent) super-Poissonian 
noise was theoretically predicted for a quantum dot sys- 
tem in the (strong inelastic) cotunneling regime analo- 
gous to ours by Sukhorukov et al44 The divergence of 
the current noise is explained as a slow switching between 
two or more current channels carrying different currents. 
We expect that the different current channels are formed 
from different resonant quantum states connecting the 
left and right dots in the cotunneling regime. Due to the 
small damping rate the switching between those channels 
is slow giving rise to the highly super-Poissonian noise. 

We also observed a quasi-divergent Fano factor (up to 
F « 600) around the shuttling instability transition point 
in the quasiclassical limit of the original one-dot shuttle 
setup ^ The explanation of the divergence is again the 
same, i.e. the slow switching between different current 
channels. Contrary to the present case the two channels 
of the one-dot setup are both given by real sequential 
tunneling processes via the dot differing just by the state 
of the oscillator (equilibrated vs. shuttling). The switch- 
ing rate between the channels can be calculated semi- 
analytically thus quantitatively confirming the proposed 
mechanismi^ In the three-dot case the semi-analytic the- 
ory would be much more complicated and we do not at- 
tempt it. A similar mechanism for the quasi-divergent 
Fano factor in a single-electron-transistor NEMS was also 
proposed recently by Blanter et al,^-^ 

Further insight to the details of the microscopic 
transport mechanism can be gained by studying the 
Wigner functions which describe the oscillator phase 
space quasiprobability distributions. We define Wigner 
functions of the unoccupied (Wuu), occupied {Wcc) cen- 
tral dot and their sum (Wtot), respectively: 

Wuu{X,P)= f ^e^""' 

X - I I {plT + pTl + P«K ) 1^ + |) ' 
Wcc{X,P) = + 
Wtot {X, P) = Wcc {X, P) + Wuu {X, P) . 

(45) 

The behavior in the close vicinity of et ~ 4?ia;o character- 
ized by a weak damping dependence of the mean current 
and the Fano factor (of the order of 1) seen in Fig. 0] 
is characteristic of shuttling. It is confirmed directly by 
the phase space plots in Fig. |S1 where the crossover from 
the predominantly shuttling transport at Sb = 3.96hujQ 
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to the cotunneling regime at Sb = 4.12fi(jJo is shown. The 
shutthng is evidenced by the asymmetric Wigner distri- 
butions of the occupied or empty central dot Wcc, Wuu, 
respectively (first column). The cotunneling manifests 
itself by the striking absence of any occupation of the 
central dot (last column) which proves the virtual nature 
of the transport in that case. 



C. Weak inter-dot coupling: sequential tunneling 
assisted by equilibrated oscillator 

Here we examine the behavior of the system in the 
weak tunneling regime, i.e. when the hopping elements 
tL^x), tji(x) coupling the adjacent dots in the array are 
small and the time scale between tunneling events is cor- 
respondingly the largest in the problem. In this limit the 
phonon subsystem gets equilibrated between the consecu- 
tive tunneling events and the distribution of the oscillator 
and bath may be taken at equilibrium corresponding to 
the appropriate electronic state. We can then solve the 
GME (O using perturbation theory keeping only the low- 
est order terms in the bare hopping parameter Vq which 
turns out to be equivalent to the P(£')-theoryiS The co- 
herence of the electron transfer process from the left to 
the right dot is broken during the transfer by the long 
enough interaction with the phonon subsystem acting as 
equilibrated thermal bath and, therefore, the resulting 
picture is just sequential tunneling (ST) via the central 
dot, at least in the device bias range where the above 
assumptions hold. We defer a more detailed discussion 
until the end of this subsection where the assumptions 
will be reexamined and their validity clarified in view of 
the obtained results. 
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FIG. 6: The four states (device empty), L (left dot occu- 
pied), C (center dot occupied), R (right dot occupied) and 
the transition rates as described by the Markov process given 
by the transition matrix 14611 . 



When we carry out the approximate solution of Eq. ^ 
in the lowest order in Vo a-s described in Appendix ^ 
we obtain the rate equation l|B6p describing a classi- 
cal Markov process of the sequential electron transfer 
between the 4 states which is depicted in Fig. El Af- 
ter introducing the vector of occupation probabilities 
p = [PqtPltPciPr\^ the equation can be rewritten in 



the matrix form p = Mp with the transition matrix 
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-{Tlc + Trc) 
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(46) 

The rates entering the matrix are calculated as func- 
tions of the model parameters from the microscopic 
P(£^)-theory and the results are given in Appendix 



Eqs. ljBT2|l . jBTSj, ljBT8|l . The stationary state p"*^* sat- 
isfying Mp"*"* = is found to be 



N 



I ^CL^RC 

rflcr + ric(rci?, + r) 
. rcL(rcfl + r) 

V ^CL^RC 



(47) 



with the normalization constant N = \V rc^ -\- 

^Lc{^cR + r) + VcL^^cR + 2rflc + r))"\ 

To calculate the mean current and, in particular, the 
current noise one can proceed following two possible 
equivalent ways which parallel in close analogy the two 
methods used in subsections [III Bl and IIII CI In the first 
method found in Refs. lsstlSa . 64 one defines an effective 
operator for the current running between, e.g. L and C 
by 




































(48) 



and together with the definition of the trace of a vector 
V as the sum of its elements, i.e. Tr v = ^ 
steady state current / reads 



J Vj , the mean 



/ = (Icl) = Tr(IcLP'*"') = iVrrciT 



R.C- 



(49) 



Using the current operator we consider the current- 
current correlation function 



CcL.CLir) = (Icl(t)Icl(O)) - (Icl)', 



(50) 



with the current-current correlator given by Hershfield et 
as 



al^ 



(Icl(t)Icl(O)) = 0(T)Tr(lcLT(r)IcLP^*'^') 
+ (?(-T)Tr(lcLT(-T)IcLP' 
+ eJ(r)Tr|lcLp''*'^*| 



stat\ 



(51) 



with the time propagator T(r) = exp(Mr) and Tr|v| = 
V'i\- This fully classical formula bears some formal 
resemblance to the quantum case (|14() but there is an im- 
portant difference in the presence of the (5-function term 
in (|51|l . While the first two terms of H51I) correspond 
to correlations between different tunneling events, the 
third term describes the self-correlation of a single tun- 
neling event within the classical description. The self- 
correlation term cannot be derived within the rate equa- 
tion formalism and was inserted by hand into the noise 
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formula of Ref. 56 based on the results of the previous 
more microscopic study,^^ Following the same line of ar- 
guments as in Sec. IIII 61 we get the following expression 
for the Fano factor F = ^ 

el 



F 



-2Tr(IcLQM-iQ IclP''*"*) + eTrllciP^ 



e(l 



CL) 



(52) 

with the projector Q = 1 - p"''^* ® [1, 1, 1, 1], = Q. 
Therefore, the Fano factor is determined by the pseu- 
doinverse of the transition matrix QM^^Q in analogy 
with the quantum-mechanical case. 

Exactly the same formula for the Fano factor can be 
obtained by employing the full counting statistics ap- 
proach analogous to the calculations in Sec. IIII CI ap- 
plied to the classical rate equation. To calculate the 
noise one has to introduce the counting variable n de- 
scribing the number of electrons that tunneled across a 
chosen junction, e.g. the LC-junction between the left 
and the central dot. Since in the present setup elec- 
trons can tunnel in the backwards direction, i.e. from 
the central dot to the left dot (see Fig.O, n can become 
negative as well. This technical detail slightly modifies 
the derivation which, however, closely follows the pre- 
vious lines. We start with Eq. H22|l where the proba- 
bility that n electrons tunneled across the LC junction 
(positive n corresponds to the left-to-center direction) 
P„(t) = Pg(")(t) + '(t) + P^'\t) + P'j^\t) is deter- 
mined by the n-resolved form of the rate equation 



p(^) F P 

^0 — ^0 
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+ rRc)P, 
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"crPr 
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(53) 



which is an intuitive generalization of the original rate 
equation ljB6|) obtained by including the transferred 
charge statistics across the LC junction, see Fig. El Per- 
forming the calculation of the noise from (|22|l in the spirit 
of Sec. IIII CI we come to the formula (|52f) again. We want 
to stress that using this second way of derivation gives 
us the entire formula with the self-correlation term and 
even the definition of the current operator H48I) appear- 
ing naturally in the course of the derivation. In this sense 
the intuitive generalization of the rate equation incorpo- 
rating the transferred charge resolution yields the full 
microscopic description of the whole process (contrary to 
the bare rate equation) and no heuristic arguments are 
necessary to get the self-correlation term. 

For the process determined by the rate matrix H4()l) the 
Fano factor can be rather easily evaluated analyticallyiS^ 
The resulting expression is, however, complicated and 
will not be given here. In the limit when F 3> 
^cLt^lc^rc^CR only the left or the central dot are 
occupied since the right dot and unoccupied state are 
immediately emptied in favor of the left dot. Due to the 




FIG. 7: The Fano factor in the two-state sequential tunneling 
limit (zero temperature, large F). The thick line is the com- 
puted Fano factor while the thin lines with circles are given by 
the formula np-|-(l — nc)'^ where nc are the occupation of the 
central dot. The collapse of the two curves marks the sequen- 
tial tunneling region. The val ues of the other parameters are 
Vo — O.lftciJo, 2:0 ~ 7 = O.Iljo, F — O.lciJo, T = 0. 



zero occupation of the right dot, the rate Fc^ despite its 
non-zero value drops out from the expressions for the sta- 
tionary probability distribution, mean current, and Fano 
factor. If, moreover, the temperature is zero we expect 
the rate Vlc to vanish (for T = only the positive de- 
vice bias range > is interesting from the ST point of 
view) and the stationary probability, mean current and 
Fano factor assume the well-known form for a two-state 



1 Z'r" 

„stat _ I RC 

^CL^RC 



F, 



CL + 



RC 
^RC 



CL -I- Tflc)^ 



(54a) 

(54b) 
(54c) 



As a consequence of these relations the Fano factor can 
be expressed in the limit F — > cx), T = in terms of, 
e.g., the stationary occupation nc — P^^^ of the central 
dot as P = -I- (1 — nc)^- This is an identity relating 
the Fano factor and the central dot occupation in the 
ST regime regardless of the particular values of the rates 
provided that the above assumptions are fulfilled. 

In Fig. Owe show the Fano factor as a function of the 
device bias for small Vq, zero temperature, and three dif- 
ferent values of a calculated numerically by the method 
described in Sec. IIII El We expect the system to be in 
the two-state ST regime described above. The thick lines 
are the Fano factor calculated directly while the thin 
lines with circles show the quantity Uq -|- (1 — nc)^ with 
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FIG. 8: Comparison between the numerical rates and the 
ones calculated by the P(i5)-theory for Vb = O.lfitJo, ct = 
G.2^Jmujo/h, xo = ^^h/niLJo, 7 = O.loJo, F = O.lwo, T = 0. 
The numerical rates are calculated assuming that the two- 
state sequential tunneling picture holds which is only true for 
£b ^ l.Sfi-tJo, see Fig. |7| In that region the two results match 
almost perfectly. 



nc being the occupation of the central dot calculated 
from the numerical evaluation of the full p*''^*. We see a 
nice collapse of the two curves for roughly Sh ^ I.^Hujq 
(depending slightly on the value of a). The collapse 
marks the two-state ST region. The discrepancy around 
< Eb < l.^hhjQ is due to cotunnehng processes pre- 
vailing over the ST ones in that region of The elec- 
tromechanical coupling terms are proportional to Sb and 
Vo and, therefore, the heat bath consisting of the me- 
chanical degrees of freedom gets almost decoupled in the 
ST regime (small Vq) at small Eh and does not suffice to 
break the coherence of the cotunnehng processes. 

We have thus verified that the identity implied by the 
two-state ST process is satisfied by the numerical results. 
While it helped us to identify the region of ST, however, 
the mentioned identity does not depend on the values 
of the rates. In the next step we calculate the values 
of the rates Tcl,^rc from the numerical results for the 
mean current and occupation of the central dot or Fano 
factor by inverting Eqs. H54|) . plot them in Fig. |S1 and 
compare with the rates calculated semi-analytically ac- 
cording to the P(i?)-theory presented in AppendixiBl We 
see a nearly perfect match between the two approaches 
in the regime of the two-state ST. The numerical rates 
were calculated using Eqs. 154() in the whole range of £b 
and, therefore, do not represent the correct rates in the 
cotunnehng dominated regime £(, ^ 1.5/zwo- The semi- 
analytical rates also confirm the cause of the ST mecha- 
nism breakdown discussed above. The Tcl rate yielding 
the bottleneck of the ST current essentially vanishes be- 
low the ST threshold and higher order processes in Vq 
(cotunnehng) take over. 

We show a representative plot of the general ST 
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FIG. 9: Comparison between the numerics and the 
P(i5)-theory based sequential tunneling picture for Vb = 
0.05/iwo, a — 0.2-y/ma;o/ft, xo = 5\/h/mujo, 7 = O.luo, F — 
O.lojo, n = 1. Due to the non-zero temperature the two-state 
model considered in Fig.|7||H|had to be extended and there is 
a new sequential tunneling region also for a negative bias. We 
observe a nearly perfect match between the two approaches 
for \eb\ ^ 1.5hujo. The behavior around = is clearly gov- 
erned by the physics of the static array since the oscillator is 
largely decoupled from the electronic degrees of freedom. 



results without the assumptions T = 0, F ^ 
^CLi^LCi^RC^CR in Fig. O The comparison between 
the numerically calculated and semi-analytical results is 
shown for both the mean current (log scale) and the Fano 
factor. Since the temperature is non-zero there is a new 
ST region for a negative bias. We see a good match be- 
tween the two approaches in the bias range |e6| ^ l.bhujQ. 
The fine structure around Sb being an even multiple of 
the oscillator frequency is given by the interplay between 
the values of different tunneling rates in those regions 
similar to the switching of the relative magnitude of Tcl 
and Tuc in Fig- HI The behavior around £{, = is clearly 
given by the physics of the static array also shown in 
the figure so that there are only small regions around 
Sb — dzhujo which are not covered either by the ST or 
static picture. To summarize, we have shown that the 
electronic transport through the array in the small Vq 
limit can be successfully described (in the device bias 
range \eb\ ^ 1.5huJo) by the ST theory with the trans- 
fer rates determined semi-analytically by the microscopic 
P(£;)-theory. 
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CONCLUSIONS 



APPENDIX A: ARNOLDI ITERATION 



We have developed theoretical techniques to evaluate 
the zero-frequency current noise in nanoelectromechani- 
cal systems. Two parallel lines have been developed: (i) 
Quantum regression theorem (QRT) and (ii) Full count- 
ing statistics (FCS). QRT has the advantage of being 
applicable to any correlation functions involving exclu- 
sively system operators, while FCS gives perhaps a more 
direct access to the current noise, but, on the other hand, 
other correlation functions cannot directly be accessed 
with it. We have demonstrated the equivalence of the 
two approaches for the model considered in this work, 
but we emphasize that the equivalence is critically de- 
pendent on whether charge conserving approximations 
are used. The three-dot model considered in this pa- 
per has a rich phenomenology allowing one to study the 
effect of the internal coherence of the electronic states, 
and by tuning the system parameters we can study the 
transition from a co-tunneling dominated regime to a se- 
quential tunneling regime. The generalized master equa- 
tions studied in this paper involve large matrices, and 
we have discussed in detail the numerical schemes that 
are needed in their solution. In certain limiting cases ap- 
proximate (semi) analytic theories can be developed, and 
we have found an excellent agreement with the full nu- 
merics. We have interpreted the computed current and 
noise curves in terms of physical concepts, and gained 
an understanding of when one can expect either sub- or 
super-Poissonian behavior. We believe that a successful 
interpretation of numerical results requires a simultane- 
ous analysis of several quantities such as mean current, 
Fano factor and Wigner distributions. 

There are several lines along which the present ap- 
proach can be continued. An interesting and impor- 
tant issue concerns the finite-frequency noise, and we are 
presently examining extensions of our theory in that di- 
rection. Spin-degree of freedom has been neglected in our 
calculations, and more work in that direction is called for. 
We have pointed out certain restrictions in the derivation 
of the generalized master equations, and one should look 
carefully at effects of (i) a more realistic confining po- 
tential, (ii) the interplay of the two different baths, and 
(iii) issues related to charge conservation. We also expect 
to get inspiration from experimental studies of quantum 
shuttles, which we hope are soon realized^ 



Acknowledgments 

We thank A. Armour, T. Brandes, T. Eirola, K. Flens- 
berg, G. Kiefilich, A. Pomyalov, and A. Wacker for stim- 
ulating discussions and comments. We are very grateful 
to A. Donarini for sharing numerical codes with us and il- 
luminating discussions. This work was supported by the 
Oticon Foundation (C.F.) and the grant 202/01/D099 of 
the Czech grant agency (T.N.) which we gratefully ac- 
knowledge. 



The key concept of the Arnoldi iteration is the 
construction of the Krylov subspace /Cj(£,xo) = 
span(xo, >Cxo, >C^xo, . . . ,£-'~^Xo) for a chosen initial su- 
pervector Xq and successively the computation of an or- 
thonormal basis Qj = [qi,...,qj] in it by the Gram- 
Schmidt orthogonalization. In the orthogonalization pro- 
cess defined by the recurrence relation 
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(A2) 

: q| • £c[k, for i = 

k<j and hk+i,k = ||>Cqfc - X^jLi hi^k q.i\\2 for k = 
j. It enters the following important relation 



CQj — Qj+i ■ Hj . 



(A3) 



Before proceeding we stress the main feature of the iter- 
ative Krylov subspace methods which consists in the fact 
that the dimension of the Krylov space is considerably 
smaller than the dimension of the original space in which 
(truncated) C acts {j — 20 in our calculations compared 
to the dimension of lON^ « 20000 of the relevant part 
of the truncated superspace). The required operations 
like finding the null space or the pseudoinverse of C are 
performed approximately in the Krylov subspace only (in 
the sense specified below) which makes them very fast. 
These fast operations are then iterated in order to achieve 
the solution of the original problem. 

The first task is to calculate the stationary density 
matrix p*'*'^* from Eq. (|43|) . This means we are looking 
for the unique null vector of the superoperator C. We 
choose an arbitrary initial vector Xq (whose choice can 
be motivated by a physical guess of the stationary state 
to improve the convergence) and construct the Krylov 
subspace /Cj(£,xo) for a fixed small j. Then we look 
for a vector x = • ^, f = (fi, ■ • ■ ,Cj)'^, ll^lb = 1 in 
the subspace which minimizes the norm ||>Cx||2 in order 
to approximate the null vector. Using l|A3|l the problem 
can be reformulated as 



mm 

ikii2=i 



mm 

ll«l|2 = l 



IIQj+i-Hj-Clb 



min II H 

ll?l|2 = l 



reil2 
(A4) 
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due to the property ||Qj+i • u||2 = ||u||2 for an arbitrary 
vector u = (ui, . . . , Uj+i)'^ . 

The last step leaves us with a problem of minimiz- 
ing the norm in a j-dimensional space spanned by the 
columns of Hj which can be solved by performing the 
singular value decomposition Hj = USVt of the rect- 
angular matrix Hj. U £ C^J+^)^(i+'^) and V e C^""^ are 



unitary matrices whereas S = 



is diagonal with positive cTfe's being the eigenvalues of 

\Jiij ■ Hj sorted in the descending orderjiS i.e. (Ti > 

cr2 > . . . CTj > 0. The norm ||Hj • ^la = ||SVf • is min- 
imized by choosing for ^ the last column of V belonging 
to the smallest singular value aj of Hj, i.e. ^ = 'Vj. The 
vector X = Qj ■ Wj is then an approximate null vector of 
C. If the norm ||>Cx||2 > tol one replaces the initial guess 
xq by X and repeats the procedure. The tolerance was 
chosen as tol = lOe ||£||2 with e being the machine preci- 
sion and the norm of the Liouvillean was estimated^^ as 
||£||2 =exp(iV/log(7V)). 

To ensure the convergence of the iteration it may be 
necessary to use preconditioning, i.e. one solves >Cx ~ 
where C = Ai~^£ with a suitable operator A4~^ which 
should be as close to the pseudoinverse of £ as possible 



J 



in order to separate the zero eigenvalue from the rest of 
the spectrum of C and thus speed up the convergenceiS 
Of course, in practice one does not have a routine for 
a pseudoinverse of C and some heuristic precondition- 
ing must be used. We used as the preconditioning the 
inverse of the "Sylvester part" Cq oi C. If we write 
Cp = Ap + pA^ + BipBj then the Sylvester part is 
given by Cop — Ap + pA^ . Performing the inversion Cq^ 
amounts to solving the Sylvester equation which is a rel- 
atively fast procedure scaling with . The usage of the 
preconditioning was in our case crucial for the conver- 
gence. After the iteration reaches its end the stationary 
density matrix is obtained by imposing the unity trace 
condition to the solution, i.e. x <-> jf*-'^'^ ^ TrgysP*'*^' = 1- 

The next step is to calculate the zero-frequency cur- 
rent noise from (|42|l . 144() . The equation H44() can be 
solved iterativcly in the Krylov subspace by the general- 
ized minimum residual method (GMRes) . If xq is an ini- 
tial approximation for the solution of £x = b the Krylov 
subspace is generated by the Arnoldi iteration starting 
with the vector Tq = b — £xo and the GMRes method 
finds a vector x e Xo-l-/Cj(£, rp) that minimizes the norm 
of the residual r = b — £x. The vector x is assumed in 
the form x = xq + Qj ■ ^ and the solution that minimizes 
the norm of the residual is obtained from 



£x||2 = min ||b — /3(xo + Qj 



= min ||Qj+i • (ei/3 - Hj 



min ||ro - £Qj • Clb = min ||ro - Qj+i • • ^||2 



min ||ei/3 — Hj • ^||2, with [3 
I 



||ro||2 and ei = (1, 0, . . . , 0)' 



(A5) 



The last minimization problem is solved easily by the 
QR- decomposition of the small rectangular matrix Hj = 
UR, where U G cO+i)xi has orthonormal columns 
(UtU = I) and R e is upper triangular. If 

has full rank the solution to the minimization problem is 
obtained by solving R • ^ = /JU^ • ei. If ||b — £x||2 > tol 
the Xq , ro are replaced by x, r and the sequence of steps 
is restarted. Again, the iteration may not converge with- 
out preconditioning. We used the same preconditioning 
as in the calculation of the null vector, i.e. we solved the 
problem £g ^£x = £q ^b by the above described algo- 
rithm. In the end of the iteration we fixed the solution 
by removing any component in the direction of the null 
vector by imposing the trace condition of (|44|l . 

It has to be noted that the choice of some suitable 
preconditioning is the difhcult part of the problem and 
most probably there is no general hint how to proceed. 
Particular cases must be attempted anew based on ex- 
perience and intuition. For example, we tried to solve 
our model for some parameters with the damping kernel 
© replaced by its translationally invariant form from 
Ref. l22i The same preconditioning yielded a convergent 
iteration scheme in much restricted range of the device bi- 



ases compared to the rotating wave approximation form 
of the damping used otherwise. Also the non-zero tem- 
perature calculations converged significantly slower than 
the corresponding zero-temperature counterparts. In the 
sequential tunneling limit the non-zero temperature cal- 
culations actually failed to converge at all so that the 
data presented in Fig. |^ had to be calculated with a di- 
rect method. Fortunately, the oscillator is in that limit 
close to its equilibrium state so that we needed A'' = 15 
at maximum which made the direct calculations feasible. 

As for the implementation of the numerical algorithms 
we used MATLAB on personal computers and/or Linux 
workstations. The building blocks are handy in MAT- 
LAB including the preconditioned GMRes routine with 
restarts which solves completely the noise calculation 
part of the problem. For efficiency reasons the station- 
ary part of the code was written "from the scratch" 
within MATLAB. The memory requirements were negli- 
gible (about 10-20 MB of RAM for N up to 40) and the 
calculation for N = 40, T = for a given set of the other 
parameters lasted a few minutes on a Linux workstation, 
moderately depending on the parameters via the number 
of required iterations to reach the convergence (a factor 
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of 2-3) . As already mentioned the non-zero temperature 
calculations were much slower and could take up to an 
hour for a given set of parameters. Most of the calcu- 
lations were done for N = 25, though, since this level 
of truncation was usually sufficient as tested by compar- 
ing results with different values of N. We also checked 
occasionally that different choices of junctions for the cal- 
culation of the mean current and the noise (|42|l gave the 
same numerical results within a very high accuracy. 



APPENDIX B: MICROSCOPIC DERIVATION OF 
THE RATE EQUATION 

In this appendix we give the derivation of the rate 
equation describing the sequential tunneling regime re- 
alized in the limit of the weak inter-dot coupling Vq — > 0. 
To this end we solve the n-unresolved version of the GME 
^ using the lowest order perturbation theory in Vq. For 
small Vq the rates (proportional to Vq) are small and 
we may assume that the oscillator gets equilibrated be- 
tween individual tunneling events between the adjacent 
dots. Within these assumptions we can find a closed set 
of equations for only the occupations of the respective 
dots Pl, Pc, Pr plus the probability that the device is 
empty Po {Pl + Pc + Pb + Po ^ 1). 

These quantities defined as Pj = {I\Trosc,B^\I) (I = 
0,L,R,C) obey the following equations stemming from 

@ 

Po = -rPo + tPr 

Pl = rPo + iTTosc,Bio-LctL{x) - tL{x)acL) 

= rPo - 2lm[TTosc,Bi(7LctL{£))] 
Pc = iTToscM^CLtLix) - tL{x)aLC 
+ ^CRtRix) - tB.(x)aRc) 

= 2Im[Trosc,B(o-LciL(i))] + 2lm[Trosc,B{^RctRix))] 
Pr = -FPr + iTiosc.Bi^RctRix)) - tn{x)acR 

= -rPfi - 2lm[TTosc.B{o-RctR{x))] . 

(Bl) 

We notice explicitly that the charge (probability) con- 
servation condition Po -I- Pl -I- Pc -I- Pr = is ful- 
filled. The occupations couple to the off-diagonal ele- 
ments (TlCi (^cr. satisfying 

O-LC = -n-^^LC + ^LC-^ X + [7?osc;0-Lc] 
V 2 ZXq 

+ i{c^LLtL{x) - tL{x)&cc) + i^LRtL{x) (B2) 

<ycR = + ^'^^~2~ ^ y^osc^ ^cr\ 

+ i{crcctR{x) - tji{x)crRR) - itL{x)aLR ~ - acR ■ 

(B3) 



In the full generality, these equations would generate an 
infinite hierarchy of equations for different moments of 
the whole density matrix a. However, in the lowest order 
in Vq we can neglect the coupling to ctlr. (which is of 
higher order in Vq) and formally integrate the equations 
leading to 



<TLc(i) 



/•OO ^ 

Jo 

/"OC 

i dr 
Jo 



-i(H' +^)tj. /u- \ i(H'~4^x)T 



e-^(HL.+ ^-^)r^^^it_^-)tL{x)e 



(B4) 



and similarly for acR(t). Now, we can employ the stan- 
dard Born-Markov approximation assuming the oscilla- 
tor plus bath subsystem in local equilibrium correspond- 
ing to a given charge state, and neglecting the memory 
effects in the evolution of P/(t)'s (both assumptions are 
justified by the small Vq): 



aLL(i-r)~aosc,B(0)PL(i) 
acc(i-T)~aosc,B(^)Pc(i) 



(B5) 



<7o.c,b(A) = e-'5(^o.o-AS)/^(A), Z{X) = 

^^^), where Trosc.B means tracing 
over the oscillator and the heat bath. 



with 



Tro.c,B(e-'3(«o, 



The rate equations for the evolution of the probabilities 
are thus: 



Po = -rPo + TPr 

Pl = rPo - TclPl + TlcPc 

Pc = YclPl - {Tlc + Yrc)Pc + TcrPr 

Pr = TrcPc - [TcR + T)Pr 



(B6) 



where the F/j's, the transition rates from the state J to 
/, are given by 
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Tlc 



2Re 
2Re 
2Re 



TcB. = 2Re 



OO 

OO 



i(e-'^°-"(Tosc,B(0)tL(x)e 



-iH'r 



dre-2"e*-"Trosc3 e 



iff' 



iff' 



^ osc 2:CQ ' 



Cosc.B 



r 



(B7) 



These rates can be also obtained starting from the Fermi 
Golden Rule expression for the bath-assisted electronic 
transitions (P(£')-theory- -) bearing in mind that the 
electronic state on the right dot is broadened by ^ due 
to the coupling to the (empty) right lead which causes 

r 

the appearance of the '^'^ factors in the expressions for 

To evaluate the rates we generalize the method used 
by Braig and Flensberg24 for the a = case. The shifted 
Hamiltonian H^g^ — -^x can be eliminated by performing 
a suitable unitary transformation which is a generaliza- 
tion of the well-known polaron shift from the indepen- 
dent boson model^^ to more oscillator modes and which 
is given by the unitary operator— 



(B8) 



where I and Ij are constants to be determined so that the 
linear shift is cancelled. It was found in Ref. '23 that 

-£b 



I 



and 



£b_ 
2xn 



2 ' 'J 



(B9) 



(BIO) 



We may thus rewrite the expression for, e.g., the Tcl 
rate as 



CL 



2Re 



dre 



(Bll) 



with the expectation value (•)o = Trosc,B(» <5'osc,b(0)). 
Using the Baker-Hausdorff theorem and introducing the 
function F{T;a) = (^e'^ir)-'^Hr)p-^A-ax\ 

F*(r; a) = F{-t; a) we get 



^g^A(r)-a^(r)g-^A-a^^^ Satisfying 



(B12a) 



Similarly, for the corresponding backward rate T]^c '^6 
get 



LC 



7j; a 



(B12b) 

with the function G(r; a) = {^e~''^(r)-ax(T)^iA-ax^^^ 
The transfer rates between the central and right dot read 

^ y^2g-2a(a;o+//2) 

r 



Zn 



(--f(l-4^))^ + (i)^ 



(B12c) 



CR 



— G[uj;a)- 



(- + f(i-4^))^ + (i)^ ■ 

(B12d) 

The evaluation of the functions F{uj; a) and G(uj; a) 
follows a standard route found in textbooks (Ref. 66, Sec. 
4.3; Ref.lH Sec. 4.4), or Ref. El Ch. 20). Technically, 
the task is to evaluate a particular characteristic func- 
tion of a (multidimensional) Gaussian distribution. The 
result is again Gaussian, into which only second-order 
correlation functions enter. 

We introduce the operator A{t; a) — A{t) + 

iax{T), so that F(T;a) = (^g.^A(r:a)^~^AH0■,a)^^ g^^^ 



G{T;a) 



-iA(r; — q) iiL+(0;-Q) 



Since H' is 



quadratic in x, Xj and p, pj we may rewrite F and 

G as F{T]a) = exp (i(2i(T; a)it (O; q) - A^{T;a) - 

it2(0;a))J, G(T;a) = exp {^{2A{r; -a)AHO; -a) - 

A^(r; —a) — A^'^{0; — q:))^) and we have thus established 
that 

G(r;a) = F(t; -a) . (B13) 

The function F{t; a) can be rewritten in terms of the 
following auxihary quantity [A cx I, see Eqs. ()B8|I. ()B9|I^ 



E{T;a;l) = (i(T; a)it (O; a))^ 

= {{A{t) +iax{T)){A{0)-iax{0)))^ 



(B14) 
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We evaluate E{t; a; I) following the lines of Ref. 12^ 
where E{t; 0; I) was evaluated. The idea is to express 
the function E in terms of the retarded Green's function 

E"{T;a;l) = -J0(r)([i(r), it(o)])^, (B15) 
using the fluctuation-dissipation theorem 

E{uj; a; I) = -2Im(i?^(w; a; I)) (l + nsiuj)) (B16) 
I 



and then find E^ by solving its equation of motion in the 
Fourier space (for details of the derivation see Ref. l65|) . 
Assuming the Ohmic coupling between the oscillator and 
the heat bath, i.e. J{uj) ~ mjujf{^), we find 



E''{u;;a;l) 



2aluj / 7 

2 1+^- 

UILOq \ u 



9 9 



which coincides with the result of Ref. 1241 for a = 0. We finally arrive at the expression for the F function 

i^(r;a)=exp / — { E{uj;a;l) e'"^'' - E{uj;0;l) + E{Lu;a;0)) 
LJ-oo 27r V /J 

2T2 



exp 



(up- — ClIq)2 + 72^;' 



T 2/ . ,2 / 



.2, .,2 



(B17) 



(B18) 



The analytical structure of the F function, in particular same as in the a ~ case^^ since it only depends on the 
the power law decay for large times at zero temperature behavior of the prefactor Puj^ - + -^^4- at w ^ 0+. 
F{t) (X T-\ T ^ oo,T = 0, S = remains the 
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The influence of two "baths" on a single system may lead to 
interesting physical phenomena — see, e.g., Ref . tfig where 
the subtle interplay between a detector and bath is studied. 
Similar study in the context of NEMS has not been carried 
out yet. 

There are different conventions used in the literature. Here 
we define the spectrum as the Fourier transform of the 
correlation function without a prefactor of 2 used in, e.g. 
in Ref. 113 The Fano factor is then given as F(0) = 
A similar proof has been suggested to us by G. Kiefilich 
(private communication) . 

We have made preliminary tests of this issue in a sim- 
pler model of the dissipative double-dot system studied 
in Ref. l5ll We have found out that the zero-frequency 
noise for the inner junction between the dots differs 
from the spectra of the outer junctions which are the 
same. This difference can be understood from the rela- 
tion [A/L.ii, /3damp] / implied by that model (cf. Ref. Isil 
Eqs. (3), (4)). On the other hand, the equivalence of the 
outer junctions is a consequence of [Ml +Mr, -Cdamp] = 0. 
Remember that the block elements po/, pio with I = 
L, C, R of the density matrix are decoupled from the rest 
and are, therefore, discarded. 

Partial and preliminary results have already been reported 

in the shuttling context 4^*— 

T. Eirola, private communication 



